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Abstract. - Bistability and hysteresis of magnetohydrodynamic dipolar dynamos generated by 
turbulent convection in rotating spherical fluid shells is demonstrated. Hysteresis appears as 
a transition between two distinct regimes of dipolar dynamos with rather different properties 
including a pronounced difference in the amplitude of the axisymmetric poloidal field component 
and in the form of the differential rotation. The bistability occurs from the onset of dynamo action 
up to about 9 times the critical value of the Rayleigh number for onset of convection and over a 
wide range of values of the ordinary and the magnetic Prandtl numbers including the value unity. 
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Introduction. — An arbitrary weak magnetic field 
may either decay or be amplified by the motion of an elec- 
trically conducting fluid. The latter case is called dynamo 
effect and it is believed to be responsible for the mag- 
netic fields of cosmic objects including the sun and most 
planets [1,2]. The dynamo effect of turbulent convection 
in rotating spherical fluid shells has received much atten- 
tion in recent years [3,4] because it is the basic model for 
the generation of the magnetic fields of the Earth and 
other planets. Many numerical studies following [5, 6] 
have attempted a direct comparison with geomagnetic ob- 
servations and have been remarkably successful in repro- 
ducing some of the main properties of the geomagnetic 
field [7, 8] while others have focused on more systematic 
explorations of the computationally accessible parameter 
space e.g. [9-13]. Dynamos obtained through numerical 
simulations are usually turbulent and it is assumed that 
their time averaged properties are independent of the ini- 
tial conditions once the computations have run for a suffi- 
ciently long time. In contrast to this assumption we wish 
to demonstrate in this Letter that bistability and hys- 
teresis of convection-driven turbulent dynamos occur in 
a wide region of the parameter space. The co-cxistcnce 
of two distinct turbulent attractors is also of general in- 
terest. In contrast to the common bistability associated 
with a subcritical bifurcation from a laminar to a turbu- 



lent attractor, the co-existence of two turbulent attractors 
is a rare phenomenon in fluid dynamics [14, 15] and mag- 
netohydrodynamics. Subcritical onset of dynamo action 
is a rather common phenomenon [16-18] and such is hys- 
teresis between non-magnetic and dynamo states. But 
co-existence and hysteresis between two fully-developed 
chaotic dynamo states far above onset are, to our knowl- 
edge, reported here for the first time. That multiple tur- 
bulent states are more likely to occur in hydromagnetic 
dynamos than in non-magnetic flows is perhaps not sur- 
prising because of the additional degrees of freedom offered 
by the magnetic field. 

Formulation. — We consider a spherical fluid shell of 
thickness d rotating with a constant angular velocity f2. 
The existence of a static state is assumed with a temper- 
ature distribution Ts = To — f3d 2 r 2 /2 and a gravity field 
in the form g = —djr, where rd is the length of the posi- 
tion vector with respect to the center of the sphere. This 
form of temperature profile alludes to the possibility that 
at least a fraction of the energy available to planetary dy- 
namos is due to radiogenic heat release. In addition to d, 
we use the time d 2 /v, the temperature v /^ad and the 
magnetic flux density vi^ig) 1 ! 2 j d as scales for the dimen- 
sionless description of the problem where v denotes the 
kinematic viscosity of the fluid, n its thermal diffusivity, 
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g its density and \i its magnetic permeability. In com- 
mon with most other simulations of Earth and planetary 
dynamos [4,7], we assume the Boussinesq approximation 
implying a constant density g except in the gravity term 
where its temperature dependence is taken into account 
with a = —(dg/ AT)/ g =const. The equations of motion 
for the velocity vector u, the heat equation for the devia- 
tion from the static temperature distribution, and the 
equation of induction for the magnetic flux density B arc 
then given by 

V-tt = 0, V-B = 0, (la) 

{d t +u- V)tt + rk x u = -Vtt + 6r + V 2 tt + B ■ VB, 

(lb) 

F(a f e + «-ve) = i?r-« + v 2 e, (ic) 

S7 2 B = P m {d t B + u- VB B Vu), (Id) 

where all gradient terms in the equation of motion have 
been combined into V7r. The dimcnsionless parameters 
in our formulation are the Rayleigh number R, the Cori- 
olis number r, the Prandtl number P and the magnetic 
Prandtl number P m , 
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where A is the magnetic diffusivity. Being solenoidal vector 
fields u and B can be represented uniquely in terms of 
poloidal and toroidal components, 

u = V x (Vv x r) + Vw x r , (3a) 
B = Vx (Vh x r) + V.9 x r . (3b) 

We assume fixed temperatures at r = = 2/3 and 
r = r = 5/3 and stress- free rather than no-slip boundary 
conditions in order to approach, at least to some extent, 
the extremely low values of viscosity believed to be appro- 
priate to planetary cores [6] , 
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For the magnetic field we assume electrically insulating 
boundaries at r = r, and r = r Q such that the poloidal 
function h matches the function which describes the 
potential fields outside the fluid shell, 
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The radius ratio r 2 /r = 0.4 is slightly larger than that 
appropriate for the Earth's liquid core. This is a stan- 
dard formulation of the spherical convection-driven dy- 
namo problem [2-4] for which an extensive collection of 
results already exists [10, 11, 19,20]. The results reported 
below are not strongly model dependent. In particular, 
dynamos with stress-free and with no-slip velocity bound- 
ary conditions as well as with different modes of energy 
supply are known to have comparable energy densities 
and symmetry properties (see fig. 15 of [4]). Furthermore, 




Fig. 1: (color online) Resolution test: Time- averaged spectra 
of magnetic (circles) and kinetic (squares) energy as a function 
of the harmonic degree I for the two cases shown in fig. 2. The 
MD case is given in blue and indicated by empty symbols and 
the FD case is given in red and indicated by solid symbols. 
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Fig. 2: (color online) Time series of two different chaotic at- 
tractors are shown - a MD (left column (a,b)) and a FD dy- 
namo (right column (d,e)) both in the case R = 3.5 x 10 6 , 
t = 3 x 10 4 , P = 0.75 and P m = 1.5. The top two panels 
(a,d) show magnetic energy densities. The rest of the panels 
show kinetic energy densities in the presence of magnetic field 
(b,e) and after the magnetic field is removed (c,f). The compo- 
nent X p is shown by thick solid black line, while X t , X p , and 
Xt are shown by thin red, green and blue lines, respectively, 
and they are also indicated by squares, triangles and circles, 
respectively. X stands for either M or E. 
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Fig. 3: (color online) A MD and a FD dynamo both in the 
case R = 3.5 x 10 6 , r = 3 x 10 4 , P = 0.75 and P m = 1.5. Each 
circle of the leftmost column shows meridional lines of constant 
B v in the left half and of r sin 9dgh = const, in the right half. 
The middle column shows lines of constant B r at r — r a + 1.3. 
Each circle of the rightmost column shows meridional lines of 
constant u v in the left half and of r sin Odgv in the right half. 

aiming to retain a general physical perspective, we inten- 
tionally use a minimal number of physical parameters in- 
cluding only those of primary importance for stellar and 
planetary applications. 

Equations of motion for the scalar fields v, w, are ob- 
tained by taking r-Vx Vx and r-Vx of equation (lb) and 
equations for g and h are obtained by taking r-Vx and 
r- of equation (Id). These equations are solved numer- 
ically by a pseudo-spectral method as described in [21] 
based on expansions of all dependent variables in spherical 
harmonics for the angular dependences and in Chebychev 
polynomials for the radial dependence. Typically, calcu- 
lations are considered decently resolved when the spectral 
power of kinetic and magnetic energy drops by more than 
a factor of 100 from the spectral maximum to the cut- 
off wavelength [9]. A minimum of 41 collocation points 
in the radial direction and spherical harmonics up to the 
order 96 have been used in all cases reported here which 
provides adequate resolution as demonstrated in fig. 1 for 
two typical dynamo solutions. The dynamo solutions are 
characterized by their magnetic energy densities, 

M p = i<| V X (Vh x r) | 2 ), M t = i(| Vg x r | 2 ), 

Mv = \{\ V x (Vh x r) | 2 ), M t = \(\ Vg x r | 2 >, 

where (•) indicates the average over the fluid shell and 
h refers to the axisymmetric component of h, while h is 
defined by h = h — h. The corresponding kinetic energy 
densities E p , Et, E p and Et are defined analogously with 
v and iv replacing h and g. 

Results. — All solutions reported below are turbulent 
and typical examples of the variations in time of the en- 
ergy densities are shown in fig. 2. Apart from the obvious 
quantitative difference, an essential qualitative change in 




Fig. 4: (color online) The ratio M p /M p of dynamo solutions as 
a function of P and P/P m at R = 3.5 x 10 6 , r = 3 x 10 4 . Full 
red and empty blue circles indicate FD and MD dynamos, 
respectively. The surface is intentionally left broken near the 
transition where it is multi- valued as discussed in text. 

the balance of magnetic energy components is observed. 
The axisymmetric poloidal component M p is dominant in 
the case shown in fig. 2(a,b) while it has a relatively small 
contribution in the case of fig. 2(d,e) with the correspond- 
ing differences in the structure of the magnetic field shown 
in fig. 3. This observation is in accord with the claim 
made in [11,20,22] that, quite generally, two regimes of 
dipolar dynamos can be distinguished, namely those with 
M p < ~M p (regime MD, "Mean Dipole") and those with 
M p > M p (regime FD, "Fluctuating Dipole"). The tran- 
sition between the regimes is shown in fig. 4 as a function 
of the ordinary Prandtl number P and of the ratio Pj P m 
in the case of fixed r = 3 x 10 4 and R = 3.5 x 10 6 . The 
transition between the two distinct turbulent dynamo at- 
tractors is not gradual, but has the character of an abrupt 
jump after a critical parameter value is surpassed as dis- 
cussed below. All solutions included in fig. 4 have a pre- 
dominantly dipolar character with the ratio M p lp /M p in 
the ranges [0.95, 1] for MD dynamos and [0.45, 0.72] for 
FD dynamos. Although, the FD dynamos feature an in- 
creased contribution of higher multipoles they are still of 
geophysical relevance [20,22]. At values of the ordinary 
and magnetic Prandtl numbers somewhat lower than those 
included in fig. 4, however, hemispherical and quadrupolar 
dynamos become predominant. 

We note that the two solutions in fig. 2 and fig. 3 are ob- 
tained at identical parameter values and thus the chaotic 
attractors MD and FD co-exist in this case. Varying 
P and P m demonstrates an extended region of this co- 
existence. In fact, the transition between the MD and 
FD dynamos is achieved via hysteresis loops as illustrated 
in fig. 5 for t = 3 x 10 4 . When an MD dynamo is used 
as initial data and the Prandtl number P is decreased, 
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Fig. 5: (color online) The upper row shows the hysteresis effect in the ratio M p /M p at r = 3 x 10 4 (a) as a function of the 
Prandtl number in the case of R = 3.5 x 10 6 , P/Pm, = 0.5; (b) as a function of the ratio P/Pm in the case of R = 3.5 x 10 6 , 
P = 0.75 and (c) as a function of the Rayleigh number in the case P — 0.75, P m = 1.5. Full red and empty blue circles indicate 
FD and MD dynamos, respectively. The critical value of R for the onset of thermal convection for the cases shown in (c) is 
R c = 659145. A transition from FD to MD dynamos as P/P m is decreased in (b) is expected, but is not indicated owing to lack 
of data. The lower row shows the value of the Nusselt number at r — Vi for the same dynamo cases. Values for non-magnetic 
convection are indicated by green squares for comparison. 



solutions remain in regime MD until the critical value 
Pfd~ 0.5 is reached at which point an abrupt jump tran- 
sition to the FD regime occurs. When a FD dynamo is 
used as initial condition and P is increased the reverse 
transition occurs at the critical value Pmi>~ 2.2 as seen 
in fig. 5(a). Similarly, coexisting attractors are found as a 
function of the ratio Pj P m for P/P m < 1 and as a function 
of the Rayleigh number above the onset of dynamo action 
as seen in figs. 5(b) and (c), respectively. The solutions 
of the hysteresis loops have been computed for up to at 
least 3 magnetic diffusion times. No evidence for a tran- 
sient nature of any case has been found. In fact, in cases 
beyond the boundaries of the hysteresis loop the transi- 
tion from MD to FD regimes or vice versa takes typically 
only 0.15 magnetic diffusion times. One may notice that 
the hysteresis loop in the direction of decreasing values of 
P/P m in fig. 5(b) is incomplete. We expect a transition 
from FD to MD dynamos to occur at sufficiently low val- 
ues of P/P m , but we are presently unable to demonstrate 
it because of the high computational costs of dynamo sim- 
ulations in this region of the parameter space. 

The bistability of convection driven dynamos is the re- 
sult of two different ways in which the magnetic field 
damps the differential rotation. In the MD dynamos the 
differential rotation generated by Reynolds stresses of the 
convection columns is eliminated almost entirely by the 
strong mean magnetic field as shown in figs. 2(b) and 3. 
Only a zonal thermal wind caused by latitudinal variations 



of the temperature remains as is typical for high Prandtl 
number dynamos [11]. Because of the strong magnetic 
field the amplitude of convection is also reduced in com- 
parison with the maximum value that it reaches in the 
absence of a magnetic field. In the case of FD dynamos 
the differential rotation is still diminished, but its align- 
ment with coaxial cylindrical surfaces is preserved. The 
amplitude of convection is now more strongly fluctuating, 
but is larger on average than in the case of the MD dy- 
namos. In this way FD and MD dynamos manage to 
carry very nearly the same heat transport as is evident 
from fig. 5. This heat transport by far exceeds the time 
average of heat transport found in the absence of a mag- 
netic field. Figure 2(c,f) demonstrates that the same type 
of convection in the form of turbulent relaxation oscilla- 
tions [19] is generated when the Lorentz force is dropped 
from the equations of motion. 

While the MD dynamos are non-oscillatory with lim- 
ited fluctuations of their axisymmctric components about 
their time average, the FD dynamos usually oscillate in 
an irregular manner. A rather pure example of such os- 
cillations is shown in fig. 6. The dipolar oscillations of 
FD dynamos can be understood in terms of Parker dy- 
namo waves [20]. In general, however, the oscillations are 
much less regular and involve considerable contributions 
of quadrupolar components. 

Conclusion. — We have demonstrated the co- 
existence of two well-distinguishable chaotic attractors in 
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Fig. 6: (color online) A half period of dipolar oscillations in the 
FD case R = 3.5 x 10 6 , r = 3 x 10 4 , P = 0.75 and P m = 0.65. 
The same fields as in the first column of fig. 3 are shown. Plots 
follow clockwise from upper left with a step At = 0.0042. 

a turbulent system over a sizable region of the parame- 
ter space. The hysteresis loops occur within the ranges 
P E (0.5, 2.2) and P/P m < 1 for the value of r used here 
and from onset of dynamo action to at least 9 x R c , where 
R c is the critical value for onset of convection. This range 
is of importance since P = 1 is used in most simulations of 
convection-driven dynamos in rotating spheres [4, 7] . The 
reason for this is that in comparisons with geomagnetic 
observations difFusivities such as v and K are interpreted 
as eddy difFusivities which are supposed to represent the 
effects of the numerically unresolved scales of the turbu- 
lent velocity field. Since the diffusion of heat and momen- 
tum owing to the small scale turbulence is similar, eddy 
difFusivities are regarded as equal with the consequence 
P = 1. 

Dynamo states similar to MD and FD dynamos have 
been obtained previously. A transition between such two 
states as a function of the Rayleigh number, R, has been 
reported in [22] where no-slip boundary conditions have 
been used. A hysteresis phenomenon was not found by 
these authors. In a later paper [13] a particular case of 
bistability has been mentioned, but no further discussion 
was given. Transitions between states similar to MD and 
FD dynamos as a function of P or P/P m caused by the 
changing strength of the inertial forces have been reported 
in [11,23]. In the simulations of [13,22] and [23] no-slip 
boundary conditions for the velocity field have been em- 
ployed and it is of interest to investigate whether, as we 
expect, the hysteresis phenomenon persist as the velocity 
boundary conditions and the heating model are changed. 
Simulations of convection-driven dynamos with no-slip 
boundary conditions and driven by a basal heat flux are 
underway and will be reported in a future paper. 

While the possibility of hysteretic behavior of planetary 
dynamos in response to slow changes in their properties 
is of considerable interest, it will be difficult to obtain 
observational evidence for such a phenomenon because of 



the long time scale of magnetic diffusion. It is important, 
however, to be aware of the bistability phenomenon in the 
interpretation of numerical dynamo simulations. 
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